clear all

load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
load matrix_time1_age0_eh0_comp1.out
%load matrix_time1_eh0_comp1.out
load stock_time2_age0_eh0_comp1_anyhs0.out
load stock_time2_age0_eh0_comp1_anyhs1.out
load stock_time2_age0_eh0_comp1.out
load stock_time2_eh0_comp1_anyhs0.out
load stock_time2_eh0_comp1_anyhs1.out
load stock_time2_eh0_comp1.out
% load trans_time2_age0_eh0_comp1.out
% load trans_time2_eh0_comp1.out
load W_time2_age0_eh0_comp1.out
load W_time2_age0_eh0_comp11.out
load W_time2_eh0_comp1.out

Par1 = zeros(35,1);
Par1(1:29) = ParEst_eh0_time0_Mex;
Par1(21) = 10^ParEst_eh0_time0_Mex(21);
Par1(22:27)=1000*ParEst_eh0_time0_Mex(22:27);
Par1(30) = 500*ParEst_eh0_time1_Mex(1);
Par1(31) = 500*ParEst_eh0_time1_Mex(2);
Par1(32:35)=ParEst_eh0_time1_Mex(3:6);

% Health transitions and unhealthy proportion
phl=matrix_time1_age0_eh0_comp1(1);
plh=matrix_time1_age0_eh0_comp1(2);
prob_l=matrix_time1_age0_eh0_comp1(3);

w_f_g=W_time2_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

w_mom=W_time2_age0_eh0_comp1; % Read wage moments
meanwi_1=w_mom(1); meanwf_1=w_mom(2);   meanwi_2=w_mom(3);  meanwf_2=w_mom(4);  sdwi_1=w_mom(5); sdwf_1=w_mom(6);  sdwi_2=w_mom(7);   sdwf_2=w_mom(8);
p10wi_1=w_mom(9);  p10wf_1=w_mom(10); p10wi_2=w_mom(11); p10wf_2=w_mom(12); p25wi_1=w_mom(13); p25wf_1=w_mom(14); p25wi_2=w_mom(15); p25wf_2=w_mom(16); 
p50wi_1=w_mom(17); p50wf_1=w_mom(18); p50wi_2=w_mom(19); p50wf_2=w_mom(20); p75wi_1=w_mom(21); p75wf_1=w_mom(22); p75wi_2=w_mom(23); p75wf_2=w_mom(24); 
p90wi_1=w_mom(25); p90wf_1=w_mom(26); p90wi_2=w_mom(27); p90wf_2=w_mom(28); minwi_1=w_mom(29); minwf_1=w_mom(30); minwi_2=w_mom(31); minwf_2=w_mom(32); 
maxwi_1=w_mom(33); maxwf_1=w_mom(34); maxwi_2=w_mom(35); maxwf_2=w_mom(36);
wage_moments_data=[meanwi_1; meanwf_1; meanwi_2; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2];

stock_mom_h=stock_time2_age0_eh0_comp1_anyhs0(1:9); % Read stocks [FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h]
stock_h=stock_mom_h';

stock_mom_l=[stock_time2_age0_eh0_comp1_anyhs1(10),stock_time2_age0_eh0_comp1_anyhs1(11),...
    stock_time2_age0_eh0_comp1_anyhs1(5),stock_time2_age0_eh0_comp1_anyhs1(8),...
    stock_time2_age0_eh0_comp1_anyhs1(9)]; % Read stocks [FX_l;IX_l;NF_l;NI_l;NN_l]
stock_l=stock_mom_l'; % not in use

stock_mom_agg=stock_time2_age0_eh0_comp1_anyhs0(1:9).*(1-prob_l)+stock_time2_age0_eh0_comp1_anyhs1(1:9).*prob_l;
stock_agg=stock_mom_agg';

HHstates_data=[stock_h;stock_agg]; % [FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h;FF_agg;FI_agg;FN_agg;IF_agg;NF_agg;II_agg;IN_agg;NI_agg;NN_agg]  joint-job state

n=50;
 % solve dist of wage offers, value functions and simulate trajectory (MSM)
[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist_2(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l,iterVF]=solve_value_function_A1_1(Par1,...
    Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions,phl,plh]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,...
    wi_1,ff_2,wf_2,fi_2,wi_2,Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l);


%Contact rates
Transitions_both = {'F-N'; 'I-N'; 'N-F'; 'N-I'; 'F-I'; 'I-F'; 'p1'; 'q1';''; 'F-N'; 'I-N'; 'N-F'; 'N-I'; 'F-I'; 'I-F'; 'p2'; 'q2'};
Estimates = [Par1(1); Par1(2); Par1(3); Par1(4); Par1(5); Par1(6); Par1(13); Par1(14);NaN;Par1(7); Par1(8); Par1(9); Par1(10); Par1(11); Par1(12); Par1(15); Par1(16)];
Transition_Rates_table = table(Transitions_both, Estimates);
writetable(Transition_Rates_table, 'Transition_Rates_table.xls')


alphaf=Par1(17);
betaf=Par1(18);
alphai=Par1(19);
betai=Par1(20);
theta=Par1(21);
a_h=Par1(22);
b1_h = Par1(23); 
b2_h = Par1(24);
a_l=Par1(25);
b1_l = Par1(26); 
b2_l = Par1(27);

gamma_h = Par1(30);
gamma_l = Par1(31);
alphaf_2=Par1(32);
betaf_2=Par1(33);
alphai_2=Par1(34);
betai_2=Par1(35);

HH_avg_Inc = p50wi_1; % median wage of informal men
MWP_gamma_h = gamma_h/(exp(-theta*HH_avg_Inc));
MWP_gamma_l = gamma_l/(exp(-theta*HH_avg_Inc));
save MWP_gamma_h MWP_gamma_h
save MWP_gamma_l MWP_gamma_l
load MWP_a_h MWP_a_h
load MWP_a_l MWP_a_l


%Preference parameters
Preference = {'theta';'b1_h';'b1_l';'b2_h';'b2_l';'a_h';'a_l';'MWP a_h';'MWP a_l';'gamma_h';'gamma_l';'MWP gamma_h';'MWP gamma_l'};
Estimates = [theta; b1_h; b1_l; b2_h; b2_l; a_h; a_l; MWP_a_h; MWP_a_l; gamma_h; gamma_l; MWP_gamma_h; MWP_gamma_l];
Preferences_table = table(Preference, Estimates);
writetable(Preferences_table, 'Preferences_table.xls')

%Wage offer parameters
Wage_offer = {'alpha_formal'; 'beta_formal';'alpha_informal';'beta_informal'};
Estimates = [alphaf; betaf; alphai; betai];
Wage_offer_par_table = table(Wage_offer, Estimates);
writetable(Wage_offer_par_table, 'Wage_offer_par_table.xls')

% Generate mean LOG WAGES - F distribution 
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wf_2f=log(sum(wf_2.*ff_2));
mean_wi_2f=log(sum(wi_2.*fi_2));
mean_wage_offer_time1=[mean_wf_1f;mean_wi_1f;NaN;mean_wf_2f;mean_wi_2f];
save mean_wage_offer_time1 mean_wage_offer_time1

%Wage offer parameters
Wage_offer = {'alpha_formal';'beta_formal';'alpha_informal';'beta_informal';'';'mean_wf_1f';'mean_wi_1f';NaN;'mean_wf_2f';'mean_wi_2f'};
Estimates = [alphaf_2; betaf_2; alphai_2; betai_2;NaN;mean_wage_offer_time1];
Wage_offer_par_table = table(Wage_offer, Estimates);
writetable(Wage_offer_par_table, 'Wage_offer_par_table_2nd_stage.xls')

x=wage_moments_data;
wage_moments_data1=[x(2);x(6);NaN;NaN;...	
x(1);x(5);NaN;NaN;...
x(4);x(8);NaN;NaN;...	
x(3);x(7)];
w=wage_moments(1:8);
wage_moments1=[w(2);w(6);NaN;NaN;... 
w(1);w(5);NaN;NaN;...
w(4);w(8);NaN;NaN;... 
w(3);w(7)];
% Model fit (for paper)
Model_fit_labels = {'mff_h'; 'mfi_h'; 'mfn_h'; 'mif_h'; 'mnf_h'; 'mii_h'; 'min_h'; 'mni_h'; 'mnn_h';'mff_agg'; 'mfi_agg'; 'mfn_agg'; 'mif_agg'; 'mnf_agg'; 'mii_agg'; 'min_agg'; 'mni_agg'; 'mnn_agg';'';'';'meanwf_1'; 'sdwf_1';'';''; 'meanwi_1'; 'sdwi_1';'';'';...
    'meanwf_2'; 'sdwf_2';'';''; 'meanwi_2'; 'sdwi_2'};
ModelFit_actual=[HHstates_data;NaN;NaN;log(wage_moments_data1)];
ModelFit_model=[HHstates;NaN;NaN;log(wage_moments1)];
Model_fit_table = table(Model_fit_labels,ModelFit_actual,ModelFit_model);
writetable(Model_fit_table, 'Model_fit_table_paper.xls')
